load('hdAll.mat');
hd=hd+hd';
K=20;n=212;Iter=10;
% center=[0, 10, 20, 30, 40, 49]+1;% 51张图 mesh编号中心位置
% center=[0, 4, 7, 3, 4, 5]+1;% 51张图 mesh编号中心位置
% all=randperm(n); center=all(1:K);
center=floor(linspace(1,212,K));
C=[];

% 第一次
label=zeros(n,1);
for i=1:length(center)
    label(center(i))=i;
end

%E步
for i=1:n
    if label(i)==0
        if i<center(1)
            label(i)=1;
            continue;
        end
        if i>center(K)
            label(i)=K;
            continue;
        end
        in=center(1:end-1)<i&center(2:end)>i;
        ind=find(in==1);
        class2=[ind,ind+1];
        center2=center(class2);
        dist=hd(i,center2);
        [~,index]=min(dist);
        label(i)=class2(index);
    end
end

cost=0;
for i=1:n
    cost=cost+hd(i,center(label(i)));
end
C=[C cost];



for ii=1:Iter
    label=zeros(n,1);
    for i=1:length(center)
        label(center(i))=i;
    end
    
    %E步
    for i=1:n
        if label(i)==0
            if i<center(1)
                label(i)=1;
                continue;
            end
            if i>center(K)
                label(i)=K;
                continue;
            end
            in=center(1:end-1)<i&center(2:end)>i;
            ind=find(in==1);
            class2=[ind,ind+1];
            center2=center(class2);
            dist=hd(i,center2);
            [~,index]=min(dist);
            label(i)=class2(index);
        end
    end
    
    %混分的情况例如 1 1 1 1 2 1 1 2 2 2 1 2 2 2 2
    %M步骤
    for i=1:length(center)
        a=find(label==i);
        ai=1:length(a);
        s=[];
        for j=1:length(a)
            s(j)=sum(hd(a(j),a(ai)));
        end
        [~,ci]=min(s);
        center(i)=a(ci);
    end
    %     center=sort(center);
    
    cost=0;
    for i=1:n
        cost=cost+hd(i,center(label(i)));
    end
    C=[C cost];
    center
    
end
% C
% tabulate(label)

% T=[1:212;center(label)]';
% csvwrite(['keyframe',num2str(K),'.csv'],T)